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(N '. Abstract 



We propose a unique scheme to construct fully optimized atomic basis sets for density-functional 
calculations. The shapes of the radial functions are optimized by minimizing the spillage of the 
wave functions between the atomic orbital calculations and the converged plane wave calculations 
for dimer systems. The quality of the bases can be systematically improved by increasing the size 
of the bases within the same framework. The scheme is easy to implement and very flexible. We 



have done extensive tests of this scheme for wide variety of systems. The results show that the 
obtained atomic basis sets are very satisfactory for both accuracy and transferability. 
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I. INTRODUCTION 

The last few years have seen the development of the first-principles methods in the ap- 
plication of complex material systems containing hundreds and thousands of atoms.- - - This 
is made possible because of the so called linear scaling algorithms^ are used, by taking the 
advantages of the locality of the electronic structures.- The widely used plane wave basis is 
not suitable for the linear scaling algorithms, because of its extended nature. Instead, the 
local basis, such as atomic orbitals are the better choices. 

The advantage of atomic orbitals are two folds. First, the basis size of atomic orbitals 
is much smaller compared to other basis sets, such as plane wave, and real space mesh, 
etc. Second, the atomic orbitals are strictly localized and therefore compatible with the 
linear scaling algorithms^ for electronic calculations. However, the atomic basis sets must 
be constructed very carefully to ensure both good accuracy and transferability. Furthermore 
the quality of the basis sets should be systematically improvable in an unbiased way. 

Several schemes to construct atomic orbitals have been proposed.-^ For example, the 
atomic orbitals can be constructed by applying certain confinement potentials to the isolated 
atoms.- 1 ^ To ensure the transferability of the orbitals, one could use larger basis set, by 
using more than one radial function for each angular moment (multi-£), or by including 
higher angular moment orbitals (polar orbitals). Empirically, the multi-^ functions can 
be generated by split-valence method,- whereas the polar orbitals can be generated by 
applying electric field in addition to the confinement potential. These methods have been 
demonstrated to be effective.- Nevertheless, different level of orbitals are treated in different 
ways, and are not guaranteed to be the optimized ones.- Volker et al. proposed a way to 
systematically improve the basis that they choose the one that improves the energy most 
from a large pool of pre-selected orbitals. 2 Alternatively, Ozaki optimize shape of atomic 
orbitals adopted to different environment as part of the self-consistency cycle. 8 However, in 
this scheme, every atom must has different orbital shape even for the same element. 

In this work, we propose a unique method that allows to construct systematically improv- 
able fully optimized atomic basis sets for density-functional calculations. Unlike previous 
methods, all the orbitals (including multi-£ and polar orbitals) can be constructed in a same 
procedure. The shapes of the radial functions are optimized by minimizing the spillage of 
the wave functions between the converged plane wave calculations and those from atomic 



orbital calculations for dimer systems, and therefore no pre-assumption about the radial 
functions is needed. The scheme is easy to implement and is very flexible and efficient. We 
have done extensive tests of this scheme for wide variety of systems. The results are very 
promising, showing very satisfactory results for both accuracy and transferability. 

The rest of paper is organized as follows. In Sec. [II] we give detailed introduction of 
our scheme to construct atomic orbitals. We test the obtained orbitals by calculating the 
structural and electronic properties of wide variety of systems in Sec. 1111] including III-V and 
group IV semiconductors, and GaN, ZnO, Al, Pb and MgO etc. We summarize in Sec. IIV1 

II. METHODS 

One of the most popular ways to generate the atomic orbitals is to use atomic orbitals of 
isolated atoms in certain confinement potentials. This procedure usually gives the minima 
basis of the atom. To make the basis more complete, one has to use multi-zeta orbitals and 
polar orbitals. The multi-zeta orbitals can be obtained via a split-valence method, 7 whereas 
the polar orbitals are generated by applying constant electric fields.- Obviously, the orbitals 
are constructed in very different procedures. The quality of the orbitals are uncontrolled, 
even though they are usually good. When larger basis sets are needed for high quality 
calculations, the procedures to get the orbitals are tedious. 

We use a very different strategy to construct fully optimized atomic orbitals that are 
highly transferable. The strategy is based on minimizing the spillage of the wave functions 
between the atomic orbital calculations and the plane wave results. The spillage is a mea- 
surement of the difference between the Hilbert space spanned by a set of local basis and the 
space spanned by the "exact" wave functions of the interested states of given systems.— 1 ^ 
The spillage is defined as,— 



n=l 



where \l/ n is the plane wave calculated eigenstate and N n is the number of states of interest. 
P is a projector which spanned by all the atomic orbitals, i.e., 

£ = £ to»)<S£<^l > (2) 



where M = M (r — r M ) is the /z-th local orbital. S^ v is the overlap matrix between orbitals 
M and (j) u , i.e., 

<S^ = (0 M l<M > ( 3 ) 

where /i = {a,i, £, /,m}, in which a is the element type, i is the index of atom of each 
element type, ( is the multiplicity of the radial functions for the angular momentum /, and 
m is the magnetic quantum number. 

The spillage has been applied to analysize the quality of given atomic basis sets.—"— There 
are also some attempts to choose optimized local basis for a given system by minimizing 
the spillage value.—"— In these methods, it usually starts from certain pre-assumed orbital 
shapes with few free parameters. The spillage is then used to determine these parameters. 
In Ref. Il2|, the authors used the spillage to choose the best Slater-type orbitals or the pseudo 
atomic orbitals for a given system. However, the transferability of the atomic orbitals are 
not taken into consideration. Using similar idea, Kenny et al. took one step further and 
generated multi-C and polar orbitals.— 

Here we propose a new scheme based on the spillage formulism to generate high qual- 
ity atomic orbitals that are systematically improvable. Our method improves upon that 
of previous methods in three aspects: (i) The shape of atomic orbitals can be generated 
automatically without any pre- assumptions, (ii) The atomic basis can be systematically im- 
proved within the same framework, (hi) The transferability of the atomic bases is improved 
by carefully choosing the reference systems. 

1. The radial functions 

In our scheme, each atomic orbital is written as a radial function multiplied by a spherical 
harmonic function. The radial functions are expanded into spherical Bessel functions. The 
fi-th local orbital is M (r) = /^(r)YJ m (f), where, 

I r > r c . 

ji(qr) is the spherical Bessel function, q is chosen to satisfy ji(qr c )—0, where r c is the cut 
off radius of the radial functions. The atomic orbitals are therefore strictly zero beyond r c . 
Yi m (r) is the spherical harmonic function, in which / is the angular momentum, m is the 



magnetic quantum number. The coefficients c m are chosen to minimize the spillage of the 
reference systems via a simulated annealing method. Since the first-principles total energy 
calculation can be done once for all for the reference system, the method is very efficient. 

We use the same energy cutoff of plane wave calculation for the spherical Bessel functions. 
For pseudopotentials calculations, 15 ~ 30 spherical Bessel functions are usually good enough 
to obtain reliable atomic orbitals. 

In order to make the kinetic energy integral well defined, one needs to make the second 
derivative of the atomic orbitals continuous. This can be done by multiplying the radial 
part of the atomic orbitals by a smooth function, 13 

^(r) = l-exp[- (r ~J c 2 ut)2 ]. (5) 

In our test, we find a has little influence on the final results, we thus fix o~=0.1. 

In our scheme, we do not have to assume the shape of the atomic orbitals, therefore in 
principle we can get the fully optimized radial functions. 

2. Reference systems 

One first needs to generate the atomic basis sets for some reference systems, then use it to 
more general cases, assuming the atomic basis sets are transferable. A good reference system 
is therefore crucial for generating high quality transferable basis set. Here, we choose dimer 
systems as reference systems. The spillage of the system is defined to be the average spillage 
values of the few selected dimers.- We choose several dimers at different bonding lengths, 
covering the whole dissociate curve of the dimer. If only one dimer is used as the reference, 
it will leave a finger print into the atomic orbitals, therefore worsen their transferability. 
We find 4 or 5 dimers are enough to generate reliable local orbitals. Further increasing the 
number of dimers to the reference system as many as 20 does not significantly improve the 
results. As it will be shown in the paper, the atomic basis generated from the dimer systems 
can be used to calculate different bulk systems with high accuracy, showing remarkable 
transferability. This is very important for studying complex material system, which may 
have complex chemical environment in a single system, including defects, surfaces, alloy etc. 



3. Systematically generate atomic orbitals 

The quality of the numerical orbital basis set can be systematically improved by increasing 
the number of radial functions (multi-zeta orbitals) of given angular momentum and by 
including orbitals with higher angular momentum (polar orbitals). There are several ways 
in the literature to construct multi-zeta orbitals and polar orbitals. However, there is still 
no systematically way to generate multi-zeta orbitals and polar orbitals, in which all orbitals 
are treated in a unbiased way. 4 ' In contrast, in our scheme all orbitals can be generated with 
the same procedures. To do so, we first generate the orbitals with given angular momentum, 
which we can call level 1 orbitals. The higher level orbitals can be generated using the same 
procedure, by minimize the spillage of the remaining Hilbert space, which orthogonal to the 
space spanned by the all previously generated atomic orbitals. Taken Si DZP (i.e., double £ 
functions plus one polar orbitals basis as an example. In the first step, we generate the first 
s and p orbitals, which form a minimal basis set for Si. In the second step, we generate the 
second s and p orbitals (multi-zeta orbitals). We first orthogonalize the wave functions of 
the reference states to the atomic orbitals generated in step 1. Here we define the projector 
operator formed by level 1 orbitals as, 

^=E*ftWi ( 6 ) 

where (j)^ is the /z-th orbital of level 1 orbitals. The remaining wave functions are, 

|< 2) > = (l-P (1) )|*i 1) > (7) 

where \&\ ) is the new set of wave functions, which is orthogonal to the atomic orbitals 
generated at step 1, i.e., 

P«|^ 2 )) = 0. (8) 

We then minimize the spillage between the second s, p orbitals and the space spanned by 
|\l/n )• In step 3, we generate the d orbitals following exactly the same procedures. The 
order of the orbitlas added into the basis can be determined by choosing the orbitals that 
decrease the spillage most after the orbitals have been added to the basis. In such way, we 
can systematically generate orbitals of any multiplicity and angular momentum in a unified 
scheme. This is important if high accuracy calculations are needed. 



4- Optimize the shape 

To give an idea of how the obtained radial functions look like using the above scheme, 
we show the radial functions of the first 3 s, p, d atomic orbitals in Fig. [T] (a), (b), (c) 
respectively for the carbon atom. The orbitals are obtained by taking five carbon dimers of 
different bond lengths as reference systems. The bond length of the dimmers are chosen to 
be 1.00, 1.25, 1.50, 2.00, 2.50 A. The energy cutoff of plane wave basis calculations is set 
to 100 Ry. The radius cutoff r c is chosen to be 6 a.u.. As we can see, the radial functions 
of these atomic orbitals have many oscillations. These oscillations are unphysical and may 
lower the transferability of the atomic basis set. To get rid of the alloying oscillations of 
the orbitals, at each step after we obtain the orbitals which minimize the spillage, we add 
a procedure to optimize the shape of the radial functions as follows. We define the "kinetic 
energy" of an atomic orbital as, 

T fl (c M ) = Y,cl q q 2 /2 + K , (9) 

q 

where q satisfy ji(qr c )—0. c m are the coefficients of the spherical Bessel functions that 
normalize the atomic orbitals, i.e., (0 M |0 M )=1. k is a penalty function that 

0, S/S -KA ^ 

oo, S/S - 1 > A 

where «So is the minimal spillage value for the given orbital set, and S is the current spillage 
value for the given coefficients c m . We found A=0.002 ~ 0.005 can be sufficient to smooth 
out the atomic orbitals. The "kinetic energies" of the atomic orbitals are also minimized via 
a simulated annealing method. After the shape optimization, the final spillage values will 
be slightly larger than the minimal ones, and have little influence on the accuracy of these 
basis sets. 

The shape optimized orbitals are plotted in Figfj] (d), (e), (f), for the s, p, d orbitals 
respectively, compared to those of unoptimized orbitals. As we see the alloying oscillations 
in the original orbitals have been gotten rid of, and the shapes of the optimized orbitals are 
much smoother than the origin ones, which implies a better transferability. We calculate 
the total energies of the reference dimmers using the optimized orbitals and find that they 
offer the same accuracy for the selected dimers as the origin ones. 



III. RESULTS AND DISCUSSION 

In this section, we do intensive tests on the accuracy and transferability of the atomic 
orbital basis sets generated using the the scheme given in Sec. [Ill for wide variety of mate- 
rials, including covalent, ionic, metallic systems. The lattice constants, bulk modulus, band 
structures calculated from the atomic basis are compared to those calculated from plane 
wave basis. Especially, GaN, ZnO and Al are known to have several stable structures that 
are energetically very close to each other, which provide very good tests on the quality of 
the atomic basis sets. 

All the calculations were done using density functional methods^ 1 ^ (DFT) with local 
density approximation(LDA) in Perdew-Zunger form.— Norm conserving pseudopotentials 18 
are used in fully separate form.— Periodic boundary condition is used for solids systems, 
and the integration over Brillouin zone is replaced by sum up Monkhorst-Pack k points.— 

A. Cutoff radii of orbitals 

The range of radial function r c is one of the most important parameters of the atomic 
orbitals. Usually larger r c leads to more accurate results, but at the same time demands 
more computational time and memory. One needs to choose proper r c that balance the two 
factors. For 3D solid, the number of neighboring atoms increases very fast as rf, one has 
to use modest r c , whereas for 2D and ID systems, r c can be relatively larger. It is also 
important to balance the errors of different elements in the system.— As one shall see below, 
it is straightforward to use the spillage value as a criterion to choose a proper r c for each 
element. Using the silicon diamond structure as an example. The energy cutoff is chosen 
to be 50 Ry and the k points are chosen to be 6x6x6, which are enough to converge the 
total energies. In Figj2] the blue curve is the total energy difference between the Si DZP 
basis and the plane wave calculations, whereas the red curve is the spillage, as a function 
of r c . As we see that both the spillage value and the total energy decrease monotonically 
as r c increases and are almost on top of each other. When r c =6 a.u., the energy difference 
can be reduced to about 1.5 eV per unit cell, and decreases rapidly to about 0.25 eV per 
unit cell at r c =8 a.u.. The spillage values for these two r c are about 9xl0 -3 and 2xl0~ 3 , 
respectively. Tests on other elements show similar results, which clearly demonstrate that 



the spillage is an excellent criterion for the quality of the atomic basis set. 

In order to further show how r c affects the spillage, we show in FigfS] (a),(b) the spillage 
of different basis size changes with r c for Si and C diamond structures. As we see, if r c is 
too small, further increasing the size of the basis does not significantly improve the quality 
of basis. For example, for Si dimers, if one chose r c =7 a.u., DZP basis has spillage about 
3.5xl0~ 3 , further increasing the basis size does not lower the spillage too much. However, 
increasing r c will dramatically decrease the spillage. For r c =12 a.u., the spillage value of DZP 
basis can be as small as 8xl0 -4 . Figure [3](b) shows similar results for diamond. However, 
as it is shown in the figure, carbon DZP orbitals with r c =Q a.u. are as good as DZP orbitals 
of silicon with r c =8 a.u.. The lesson we can learn from the tests is that one should choose 
proper basis set size for a given r c . 

B. III-V and group IV semiconductors 

Semiconductors are a class of important materials. Here we test our atomic bases for 
III-V and group IV semiconductors. The energy cutoff is chosen to be 50 Ry and the k 
points are chosen as 6x6x6, unless otherwise noticed. The results are summarized in Table 
HI The r c of Ga, In, Al, As, P, Sb, Ge elements are chosen to be 9 a.u. whereas r c of 
silicon and carbon is chosen to be 8 a.u. and 6 a.u., respectively. We use plane wave basis 
results as benchmark, also listed in TableU We can see that the results from the single- 
zeta (SZ) basis (the minimal basis) have large deviations for both lattice constants and 
bulk modulus. SZ basis predict too large lattice constants than those calculated from plane 
wave basis for more than 0.1 -0.2 a.u. and underestimate the bulk modulus for more than 
10%. However, modest size double-zeta plus polarized basis (DZP) basis always offer good 
results. The largest difference occurs in Ge, where the deviation of lattice constant and bulk 
modulus between DZP and plane wave basis is 0.07 a.u. and 4 GPa, respectively. After 
increasing the number of basis to triple- zeta plus double polarized basis (TZDP, we also use 
notation 3s3p2d), we can see systematically improvement over the DZP orbitals, and are 
almost identical to those calculated from plane wave basis. The atomic orbitals generated 
from dimer reference system can also obtain such good results for solid systems, showing 
remarkable transferability. 

The band structures of silicon (diamond structure) are shown in Fig. 0] The lattice 
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constant is fixed at 10.20 a.u., the energy cutoff is 50 Ry. We compare the band structures 
calculated by atomic DZP, TZDP and 5ZQP (5s5p4d) bases in Fig. 0(a) (b)(c) respectively 
and compared with the plane wave results. The black solid lines represent the plane wave 
results, whereas the blue dotted lines are the results of atomic bases. The Fermi level is 
fixed at eV. From FigI4(a) we can see DZP basis already gives very good band structures, 
for valence bands and the low energy conduction bands, except around L point. The TZDP 
basis improves the band structures around L point, as illustrated in Fig0(b), though there 
is still small difference. If we further increase the basis size to 5ZQP, we can see the energy 
bands calculated from atomic orbitals are almost identical to those from plane wave basis. 
We also calculate the electronic states of a Silicon cluster containing 29 silicon atoms and 
38 hydrogen atoms. We calculate lowest 100 states and plot the density of states (DOS). Also 
we denote 2slp, 3s2p and 4s3p as DZP, TZDP and QZTP basis for hydrogen, respectively. 
The results are shown in Fig. [5(a), (b),(c) for the DZP, TZDP and QZTP bases respectively. 
The DOS calculated by plane wave is shown in black solid line, where as those calculated 
from atomic orbitals are shown in red dashed lines. For the DZP basis, we find that the DOS 
of the valence states are almost identical to that calculated from plane wave basis. However, 
the DOS of conduction electrons shifts to the higher energy side for about 280 meV relative 
to the plane wave result. The TZDP basis improves the DZP's results, however the DOS 
of conduction electrons still shifts a little towards the higher energy side. QZTP further 
improves the results, which are in excellent agreement with those calculated from plane 
wave basis. 

C. GaN, ZnO, Al, MgO, Pb 

Now we test our atomic bases for several important materials, including GaN, ZnO, Al 
and MgO and Pb. GaN, ZnO and Al have several stable structures. Zn has 3d electrons as 
valence electrons, whereas Al is metallic. They thus offer good examples for comprehensive 
tests on the quality and transferability of the atomic bases. 

GaN has two stable crystal structures, namely, the zinc blende structure (B3) and the 
wurtzite structure (B4). The energy difference between the two structures is very small. 
The B3 structure has only one structure parameter, i.e., the lattice constant, whereas the 
B4 structure can be described by three parameters: the lattice constants a, c and the 
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internal parameter u, which describes the relative position of the two hexagonal close-packed 
sublattices. The purpose of the test is to see if the atomic basis set can predict correctly the 
energy difference between the two structures. The energy cutoff is chosen 120 Ry and the 
k points is chosen to be 6x6x6 for B3 and 6x6x4 for B4. The calculated total energies as 
functions of volume per atom are shown in Fig. El compared to those calculated from plane 
wave basis. The plane wave calculations indicate that wurtzite structure is more stable 
than zinc blende structure, which is also predicted correctly from both DZP and TZDP 
basis sets. We can also see that the total energies are systematically improved from DZP to 
TZDP. More properties including structure parameters (a, c and u), bulk modulus and the 
total energies difference between the two structures are shown in Table [Til We can see that 
DZP basis already give very good structure parameters and bulk modulus compared with 
those calculated from plane wave basis for both B3 and B4 structures. TZDP basis further 
improves all properties. It gives much better bulk modulus than DZP basis, which reduce 
the difference from about 10 GPa to less than 1 GPa. The total energy difference between 
zinc blende structure and the wurtzite structure from plane wave calculations is about 6 
meV. Although the energy difference is very small, both DZP and TZDP basis give rather 
good value, which are 4 meV and 7 meV, respectively. This proves that for GaN, local basis 
can provide accurate results as good as plane wave basis. 

Let us see how the scheme work for systems containing 3d electrons. Zinc is a transition- 
metal element, which contains 3d electrons. We calculate four structures of ZnO, including 
rock salt structure (Bl), cesium chloride structure (B2), zinc blende structure (B3) and 
wurtzite structure (B4). The energy cutoff is chosen as 120 Ry, and we use 6x6x6 /c-meshes 
for Bl, B2 and B3 structures and 6x6x4 fc-meshes for B4 structure. We define 2s2pld for 
O and 2slp2d for Zn as DZP basis, whereas 3s3p2d for O and 3s2p3d for Zn, as TZDP basis. 
The r cut is 8 a.u. for O and 8 a.u. for Zn. We plot total energies vs volume per atom for 
these four ZnO structures in Fig. [3(b) using DZP basis, compared with the results calculated 
from plane wave basis shown in Fig. [3(a). Both bases give correct energy order for the four 
structures compared to experiments. 22 As we see the energy diagrams calculated from DZP 
basis look almost the same as those calculated from plane wave basis, except that all the 
total energies calculated from DZP basis shift up for about 0.185 eV per atom relative to 
the plane wave results. The calculations shows that the structures with decreasing energy 
order is B2, Bl, B3 and B4, and ground state structure of ZnO is the wurtzite structure. 

11 



Table [TTT] further shows the calculated structure parameters, bulk modulus and energy 
differences of ZnO using DZP and TZDP basis comparing with plane wave results. Both 
DZP and TZDP basis give accurate lattice constants. The largest difference is less than 
0.02 A compared to plane wave results. However, TZDP basis gives more accurate bulk 
modulus for all structures. The energy difference between wurtzite phase and zinc blende 
phase is also calculated accurately using different basis. The energy of structure B4 is 9 
meV per atom lower than that of the structure B3 as calculated by plane wave basis. DZP 
and TZDP give exact the same results. We also show the lattice constant and bulk modulus 
calculated by Gaussian basis. 22 We see that the errors due to pseudopotentials and other 
approximations are bigger than the errors caused by the atomic bases. 

Aluminum is a metal system. It is not obvious that atom-centered atomic basis can 
describe it well. We test four structures of Aluminum, including the simple cubic (sc) 
structure, the face-centered cubic (fee) structure, the body-centered cubic (bec) structure 
and the hexagonal-closed packed (hep) structure. The energy cutoff is fixed at 70 Ry, and 
the k points are chosen as 6x6x6. Gaussian smearing is used in all calculations. We use 
DZP (2s2pld) orbitals for Aluminum with r c =9 a.u. Figure [5(a) , (b) compare the energy 
diagrams of the Aluminum four structures calculated using plane wave basis and DZP basis. 
As we see, the atomic basis provide excellent agreement with the plane wave result for all 
structures. The fee structure is the lowest-energy structure of Aluminum predicted by both 
bases. We summarize calculated properties, including lattice constant, bulk modulus and 
energy difference between different structures of Aluminum in table HVl Surprisingly, we find 
that DZP basis can provide extremely good results compared to plane wave calculations.—. 
For example, in our calculation, the largest difference of calculated lattice constant is 0.011 
A in bec structure, whereas the largest difference of bulk modulus is 1 GPa. The DZP basis 
can also give excellent energy differences between different structures of Aluminum. These 
results are much better than previous calculations, also using DZP basis. 13 For example, in 
the lattice constants errors for fee, hcp(a/c) bec and sc are 0.044 A 0.054/0.040 A 
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0.037 A 0.013 A respectively. The bulk modulus errors for fee, hep, bec and sc is 4.9 GPa, 
8.5GPa. 5.2 GPa and 2.9 GPa, respectively, which are much larger than those obtained 
using our bases. 

We also test the quality of our atomic bases for some other materials, such as MgO, Pb, 
etc. The results are summarized in Table El The energy cutoff is 70 Ry for MgO and 50 
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Ry for Pb. The k points are 6x6x6. We all use 2s2pld for Mg, 0, Pb. We use r c =8 a.u. 
for Mg and 0, and r c =9 a.u. for Pb. We compare the results with previous calculations 
using DZP basis.— We find our basis give much better bulk modulus for Pb than previous 
calculations. 

IV. SUMMARY 

We propose a unique scheme to construct fully optimized atomic basis sets for density- 
functional calculations. The shape of the radial functions are optimized by minimizing the 
spillage of the wave functions between the atomic orbital calculations and the converged 
plane wave calculations for dimer systems. Our method improves upon that of previous 
methods in three aspects: (i) The shape of atomic orbitals can be generated automatically 
without any pre-assumptions. (ii) The atomic basis can be systematically improved within 
the same framework. (iii)The transferability of atomic orbitals bases are improved by care- 
fully choosing the reference systems. The scheme is easy to implement and very flexible. 
We have done extensive tests of this scheme for wide variety of systems, including semicon- 
ductors, ionic, covalent and metallic materials. The results show that the obtained atomic 
basis sets are very satisfactory for both accuracy and transferability. 
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TABLE I: Comparison of the calculated lattice constants a (in a.u.) and bulk modulus B (in GPa) 
using plane wave basis (PW) and atomic bases for typical III-V and group IV materials. 

a B 

Compound SZ DZP TZDP PW Expr. SZ DZP TZDP PW Expr. 



GaAs 
GaP 
GaSb 
InAs 

InP 

InSb 

AlAs 

A1P 

AlSb 

Ge 

Si 

C e 



10.67 10.50 10.49 10.48 10.68 a 69 78 

10.28 10.11 10.11 10.10 10.30 a 82 92 

11.54 11.38 11.37 11.36 11.52 a 49 59 

11.40 11.27 11.28 11.28 11.45 a 63 66 

11.07 10.94 10.94 10.93 11.09 a 78 79 

12.33 12.05 12.05 12.07 12.24 a 41 50 

10.88 10.63 10.62 10.59 10.70 a 67 76 

10.50 10.26 10.25 10.21 10.33 a 64 87 

11.83 11.58 11.57 11.54 11.59 a 48 55 

10.82 10.68 10.61 10.61 10.70 c 57 67 

10.59 10.28 10.25 10.23 10. 26 6 74 94 

6.78 6.67 6.67 6.67 6.75 b 436 470 467 466 442 fc 



77 


77 


75.57 c 


93 


93 


89 d 


58 


57 


hl d 


65 


65 


m d 


79 


80 


71 d 


49 


50 


4jd 


76 


76 


77 d 


88 


89 


86 d 


56 


56 


58 d 


71 


71 


77.20 c 


94 


94 


99 6 



a I. Vurgaftman. 
6 C. Kittel, Ref. 



24 



R. Meyer and L. R. Ram- Mohan, Ref. 



23. 



c Yu-Min Juan and Eftbirriios Kaxiras, Ref 
d Marvin L. Cohen, Ref|2J. 
e Energy cutoff 100 Ry is used. 
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TABLE II: Basis comparision for GaN Zinc blende (B3) and wurtzite (B4) structures, a, c (in A) 
are the lattice constants, and c (in A) is the internal parameter. B (in GPa) is the bulk modulus, 
whereas AE (in eV/atom) is the total energy difference between different structures. The energy 



set to zero 


for B4 structure. The data of ' 


Other calculations" 


and experim 


om the reference 27 










Compound Proper 


ies PW DZP 


TZDP Other calculations 


Experiment 


GaN(B3) 


a 


4.424 4.441 


4.422 


4.446-4.46 


4.49 




B 


207 197 


208 


184-207 


173 




AE 


0.006 0.004 


0.007 






GaN(B4) 


a 


3.130 3.142 


3.130 


3.126-3.170 


3.160-3.189 




c/a 


1.629 1.631 


1.630 


1.620-1.638 


1.621-1.626 




u 


0.377 0.377 


0.377 


0.377-0.379 






B 


207 194 


207 


190-245 


188-237 




AE 
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TABLE III: Basis comparison for ZnO rock salt (Bl), cesium chloride (B2), zinc blende (B3) and 
wurtzite (B4) structures, a, c (in A) are the lattice constants, and c (in A) is the internal parameter. 
B (in GPa) is the bulk modulus, whereas AE (in eV/atom) is the total energy difference between 
different structures. The energy difference is set to zero for B4 structure. The results from Gaussian 



orbitals and the experimental values are taken form Ref. 



22 



Compound Properties PW DZP TZDP Gaussian Experiment 

ZnO(Bl) a 4.286 4.296 4.296 4.218 4.271-4.283 

B 196 198 195 203 177-228 

AE 0.068 0.066 0.058 

ZnO(B2) a 2.653 2.659 2.658 2.614 

B 194 190 191 201 

AE 0.639 0.647 0.620 

ZnO(B3) a 4.583 4.596 4.595 4.509 4.62 

B 152 150 151 154 

AE 0.009 0.009 0.009 

ZnO(B4) a 3.256 3.266 3.264 3.205 3.248-3.250 

c 5.246 5.256 5.255 5.151 5.207-5.210 

u 0.381 0.381 0.381 0.381 0.382 

B 152 150 151 155 136-183 
AE 
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TABLE IV: Comparison of the calculated lattice constants a, c (in A), bulk modulus B (in 
GPa) and energy difference AE (in eV/atom) for Aluminum fee, bec, sc and hep strcutures, using 
different basis. The energy difference is set to zero for fee structure. The results of "Other DZP" 
and "Other PW" are taken from Ref. 13. 



Compound Properties 


PW 


DZP 


Other DZP Other PW 


Al(fcc) 


a 


3.964 


3.974 


4.011 


3.967 




B 


80 


79 


76.2 


81.1 




AE 














Al(bcc) 


a 


3.175 


3.186 


3.212 


3.175 




B 


75 


75 


67.1 


72.3 




AE 


0.1167 0.1157 


0.128 


0.121 


Al(sc) 


a 


2.678 


2.678 


2.688 


2.675 




B 


54 


54 


57.9 


60.8 




AE 


0.5819 0.6018 


0.479 


0.483 


Al(hcp) 


a 


2.806 


2.805 


2.846 


2.793 




c 


1.638 


1.639 


1.63 


1.67 




B 


82 


83 


71.6 


80.1 




AE 


0.0836 0.0821 


0.032 


0.038 
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TABLE V: Comparison of the calculated lattice constants a (A) and bulk modulus B (GPa) using 
atomic orbitals to those from plane wave calculations for MgO and Pb. The data of "Other DZP", 
"Other PW" and experimental values are taken from Ref. lid . 

Compound Properties PW DZP Other DZP Other PW Experiments 



MgO(Bl) 


a 


4.233 4.238 


4.11 


4.10 


4.21 




B 


169 169 


167 


168 


152 


Pb(fcc) 


a 


4.873 4.875 


4.88 


4.88 


4.95 




B 


56 55 


64 


54 


43 



19 



Minimize spillage only Minimize kinetic energy 




1234560 123456 
Radius (a.u.) Radius (a.u.) 



FIG. 1: (Color online) (a), (b), (c) The radial functions of unoptimized carbon s, p, d orbitals 
respectively, (d), (e), (f) The radial functions of optimized carbon s, p, d orbitals. 
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FIG. 2: (Color online) The total energy difference (blue line) and average spillage value of five 
dimers (red line) as functions of orbital radius cutoff r c for Si DZP orbitals. 
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FIG. 3: Convergence of spillage value as functions of orbital radius cutoff r c and basis size for (a) 
Si and (b) Carbon dimers. 
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FIG. 4: (Color online) Comparision of the band structures of Si calculated by different atomic 
bases (blue dotted lines) and plane wave basis (solid balck lines). 
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FIG. 5: (Color online) Comparison of the density of states (DOS) of S129H38 cluster calculated by 
different atomic bases and plane wave basis. 
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FIG. 6: (Color online) Comparison of the total energies of GaN zinc blende structure (B3) and 
the wurtzite structure (B4) as functions of volume per atom using different atomic bases to those 
using plane wave basis. 
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FIG. 7: (Color online) Comparison of the total energies of ZnO rock salt structure (Bl), cesium 
chloride structure (B2), zinc blende structure (B3) and wurtzite structure (B4) as functions volume 
per atom (a) using plane wave basis to (b) those using atomic DZP basis. 
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FIG. 8: (Color online) Comparison of the total energies of Al sc, fee, bcc, hep structures as functions 
of volume per atom (a) using plane wave basis to (b) those using atomic DZP basis. 
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